///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file restricted_pairing_function.cpp
///
/// \brief   implementation for the restricted pairing function class
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#include<src/pcps.h>
#include<src/input.h>
#include<src/coeff_obj.h>
#include<src/mantissa_exponent.h>
#include<src/wfn_data.h>
#include<src/hamiltonian.h>
#include<src/command.h>
#include<src/restricted_pairing_function.h>
#include<src/matrix_math.h>
#include<src/blas_lapack.h>
#include<src/timing.h>

// ///////////////////////////////////////////////////////////////////////////////////////////////////
// /// \brief returns the value of the coefficient for the currently initialized configuration
// ///        with noise removed
// ///
// /// \param[in]      config   the current lattice and noise configuration
// ///
// ///////////////////////////////////////////////////////////////////////////////////////////////////
// template<class S> pcps::MantExp<S> pcps::RestPairingFunc<S>::noiseless_value(const pcps::Configuration<S> & config) {
//   
//   // compute noiseless value via matrix determinant lemma
//   const S x11 = -pcps::xdot(_nap, &config.noise_vecs_v().at(0).at(0), 1, &_Iu1.at(0), 1);
//   const S x12 = -pcps::xdot(_nap, &config.noise_vecs_v().at(0).at(0), 1, &_Iu2.at(0), 1);
//   const S x21 = -pcps::xdot(_nap, &config.noise_vecs_v().at(1).at(0), 1, &_Iu1.at(0), 1);
//   const S x22 = -pcps::xdot(_nap, &config.noise_vecs_v().at(1).at(0), 1, &_Iu2.at(0), 1);
//   return this->_curr_val * ( ( pcps::unity<S>() + x11 ) * ( pcps::unity<S>() + x22 ) - x21 * x12 );
// 
// }
// 
// template pcps::MantExp<double> pcps::RestPairingFunc<double>::noiseless_value(const pcps::Configuration<double> &);
// template pcps::MantExp<std::complex<double> >
//          pcps::RestPairingFunc<std::complex<double> >::noiseless_value(const pcps::Configuration<std::complex<double> > &);

// ///////////////////////////////////////////////////////////////////////////////////////////////////
// /// \brief Tests whether a rank two matrix determinant lemma correctly removes the noise
// ///
// /// \param[in]      config   the lattice configuration
// ///
// ///////////////////////////////////////////////////////////////////////////////////////////////////
// template<class S> void pcps::RestPairingFunc<S>::test_noise_removal(const pcps::Configuration<S> & config) {
//   
//   // compute noiseless value via matrix determinant lemma
//   const S noiseless = this->noiseless_value(config).get();
// 
//   // compute noiseless value directly
//   pcps::Configuration<S> test_config(config);
//   for (int i = 0; i < _nap; i++) {
//     test_config.noise_vecs_u().at(0).at(i) = pcps::zero<S>();
//     test_config.noise_vecs_u().at(1).at(i) = pcps::zero<S>();
//     test_config.noise_vecs_v().at(0).at(i) = pcps::zero<S>();
//     test_config.noise_vecs_v().at(1).at(i) = pcps::zero<S>();
//   }
//   pcps::RestPairingFunc<S> test_func(*this);
//   test_func.check_initialize(test_config, false);
// 
//   // print a comparison of the two values
//   std::cout << boost::format(" determinant lemma = %16s") % pcps::format_number("%16.12f", noiseless) << std::endl;
//   std::cout << boost::format("direct calculation = %16s")
//                % pcps::format_number("%16.12f", test_func._curr_val.get()) << std::endl;
//   std::cout << boost::format("             error = %16s")
//                % pcps::format_number("%16.12f", noiseless - test_func._curr_val.get()) << std::endl;
//   std::cout << std::endl;
// 
// }

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  prints the descriptor.
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void pcps::RPF_Descriptor<S>::printout(const std::string & name) const {

      std::stringstream a_sites_ss;
      std::stringstream b_sites_ss;
      std::stringstream  a_perm_ss;
      std::stringstream  b_perm_ss;
      std::stringstream  a_ppos_ss;
      std::stringstream  b_ppos_ss;
      std::stringstream move_ss;
      for (int i = 0; i <        _a_sites.size(); i++) a_sites_ss << " " <<        _a_sites.at(i);
      for (int i = 0; i <        _b_sites.size(); i++) b_sites_ss << " " <<        _b_sites.at(i);
      for (int i = 0; i < _a_ordered_perm.size(); i++)  a_perm_ss << " " << _a_ordered_perm.at(i);
      for (int i = 0; i < _b_ordered_perm.size(); i++)  b_perm_ss << " " << _b_ordered_perm.at(i);
      for (int i = 0; i <     _a_perm_pos.size(); i++)  a_ppos_ss << " " <<     _a_perm_pos.at(i);
      for (int i = 0; i <     _b_perm_pos.size(); i++)  b_ppos_ss << " " <<     _b_perm_pos.at(i);
      move_ss << _old_occ_site << " --> " << _new_occ_site;
      std::cout << boost::format("%6s  %16s  %16s  %16s  %16s  %16s  %16s  %16s")
                   % name
                   % a_sites_ss.str()
                   % b_sites_ss.str()
                   %  a_perm_ss.str()
                   %  b_perm_ss.str()
                   %  a_ppos_ss.str()
                   %  b_ppos_ss.str()
                   %    move_ss.str()
               << std::endl;

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  Constructs the descriptor.
///
/// \param[in]      userinp     the user's input options
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> pcps::RPF_Descriptor<S>::RPF_Descriptor(const pcps::Input & userinp)
  : _nas(userinp.nsites()/2),
    _nap(userinp.nparticles()/2),
    _new_occ_site(-1),
    _old_occ_site(-1),
    _a_sites(userinp.nparticles()/2, 0),
    _b_sites(userinp.nparticles()/2, 0),
    _a_ordered_perm(userinp.nparticles()/2, 0),
    _b_ordered_perm(userinp.nparticles()/2, 0),
    _a_perm_pos(userinp.nparticles()/2, 0),
    _b_perm_pos(userinp.nparticles()/2, 0)
{}

template pcps::RPF_Descriptor< double               >::RPF_Descriptor(const pcps::Input &);
template pcps::RPF_Descriptor< std::complex<double> >::RPF_Descriptor(const pcps::Input &);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief Initializes the permutations of the rows or columns that would put them in ascending
///        order based on their sites
///
/// \param[in]      sites       the sites corresponding to the rows or columns
/// \param[out]     ordered     on exit, the permutation that would put the rows/columns in order
/// \param[out]     perm_pos    on exit, the position of each row/column in the permutation
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S> void pcps::RPF_Descriptor<S>::initialize_permutation(const std::vector<int> & sites,
                                                                       std::vector<int> & ordered,
                                                                       std::vector<int> & perm_pos)
{

  assert( sites.size() == _nap );
  assert( ordered.size() == _nap );
  assert( perm_pos.size() == _nap );

  int lower_bound = -1;
  for (int i = 0; i < _nap; i++) {
    int min_val = _nas;
    int min_ind = -1;
    for (int j = 0; j < _nap; j++) {
      if ( sites[j] < min_val && sites[j] > lower_bound ) {
        min_val = sites[j];
        min_ind = j;
      }
    }
    assert( min_ind >= 0 );
    assert( min_ind < _nap );
    lower_bound = min_val;
    ordered[i] = min_ind;
    perm_pos[min_ind] = i;
  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  Prepares the various arrays that track the pairing function's occupancy information
///
/// \param[in]      config   the lattice configuration
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S> void pcps::RPF_Descriptor<S>::prepare_occ_info(const pcps::Configuration<S> & config) {

  // populate the alpha and beta site lists
  int a_count = 0;
  int b_count = 0;
  for (int i = 0; i < 2*_nas; i+=2) {
    if (config[i] == 1)
      _a_sites.at(a_count++) = i/2;
    if (config[i+1] == 1)
      _b_sites.at(b_count++) = i/2;
  }
  assert(a_count == _nap);
  assert(b_count == _nap);

  // initialize the permutation arrays
  this->initialize_permutation(_a_sites, _a_ordered_perm, _a_perm_pos);
  this->initialize_permutation(_b_sites, _b_ordered_perm, _b_perm_pos);

  // initialize so there is no move
  _new_occ_site = 2 * _a_sites.at(0) + 0;
  _old_occ_site = 2 * _a_sites.at(0) + 0;

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief Updates the descriptor for a configuration in which one electron has moved.
///
/// \param[in]      new_ind     the new position of the electron
/// \param[in]      old_ind     the old position of the electron
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S> void pcps::RPF_Descriptor<S>::update(int new_ind, int old_ind) {

  // get whether the move was alpha or beta
  const int b = old_ind % 2;
  assert( new_ind % 2 == old_ind % 2 );

  // get the row or column index of the destroyed electron
  const int rc_ind = this->get_row_col_index(old_ind);

  // update the occupancy list
  std::vector<int> & occ_sites = ( b == 0 ? _a_sites : _b_sites );
  assert( occ_sites.at(rc_ind) == old_ind / 2 );
  occ_sites.at(rc_ind) = new_ind / 2;

  // update the permutation that keeps track of the ordering of the rows' or columns' sites
  if ( b == 0 )
    this->initialize_permutation(_a_sites, _a_ordered_perm, _a_perm_pos);
  else
    this->initialize_permutation(_b_sites, _b_ordered_perm, _b_perm_pos);

  // initialize so there is no move
  _new_occ_site = 2 * _a_sites.at(0) + 0;
  _old_occ_site = 2 * _a_sites.at(0) + 0;

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   builds the occupied pairing matrix based on the current occupancy information
///
/// \param[in]      pm      the full pairing matrix
/// \param[out]     opm     on exit, the occupied pairing matrix
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S> void pcps::RPF_Descriptor<S>::build_occupied_matrix(const std::vector<S> & pm, std::vector<S> & opm) const {

  // ensure the matrix is correctly sized
  if (opm.size() != _nap*_nap)
    throw pcps::Exception("matrix size (%i) is wrong in pcps::RPF_Descriptor<S>::build_occupied_matrix") % opm.size();

  //for (int i = 0; i < _nap; i++)
  //  std::cout << _a_sites.at(i) << " ";
  //std::cout << std::endl;
  //std::cout << std::endl;

  //for (int i = 0; i < _nap; i++)
  //  std::cout << _b_sites.at(i) << " ";
  //std::cout << std::endl;
  //std::cout << std::endl;

  // build the bare occupied pairing matrix
  for (int i = 0; i < _nap; i++)
  for (int j = 0; j < _nap; j++) {
    assert( i*_nap+j >= 0 && i*_nap+j < opm.size() );
    assert( _a_sites[i] * _nas + _b_sites[j] >= 0 && _a_sites[i] * _nas + _b_sites[j] < pm.size() );
    opm[i*_nap+j] = pm[ _a_sites[i] * _nas + _b_sites[j] ];
  }

  //// add noise to the pairing matrix
  //const std::vector<S> & u1 = config.noise_vecs_u().at(0);
  //const std::vector<S> & u2 = config.noise_vecs_u().at(1);
  //const std::vector<S> & v1 = config.noise_vecs_v().at(0);
  //const std::vector<S> & v2 = config.noise_vecs_v().at(1);
  //for (int i = 0; i < _nap; i++)
  //for (int j = 0; j < _nap; j++) {
  //  assert( i*_nap+j >= 0 && i*_nap+j < opm.size() );
  //  assert( i >= 0 && i < u1.size() );
  //  assert( i >= 0 && i < u2.size() );
  //  assert( j >= 0 && j < v1.size() );
  //  assert( j >= 0 && j < v2.size() );
  //  opm[i*_nap+j] += u1[i] * v1[j] + u2[i] * v2[j];
  //}

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  gets the row or column index of the given site.
///
/// \param[in]      site      the site whose row or column index we wish to find
///
/// \return the row or column index of the site
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> int pcps::RPF_Descriptor<S>::get_row_col_index(int site) const {

  // get the appropriate vector of occupied sites
  const std::vector<int> & occ_sites = ( site % 2 == 0 ? _a_sites : _b_sites );

  // change from all site indexing to alpha or beta only indexing
  site = site/2;

  // look for the site
  int retval = -1;
  for (int i = 0; i < _nap; i++) {
    if ( occ_sites.at(i) == site ) {
      if (retval >= 0)
        throw pcps::Exception("site appeared more than once in pcps::RPF_Descriptor<S>::get_row_col_index");
      retval = i;
    }
  }

  // check that the site was found
  if (retval < 0)
    throw pcps::Exception("site not found in pcps::RPF_Descriptor<S>::get_row_col_index");

  // return the site
  return retval;

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief computes the difference between the new and old row or column for a one electron
///        change to the fock space configuration.
///
/// \param[in]      pm          Size _nas*_nas.  The full pairing matrix.
/// \param[out]     vec         Size _nap.  On exit, the row or column difference.
/// \param[in]      new_ind     the index of the newly occupied site (defaults to _new_occ_site)
/// \param[in]      old_ind     the index of the newly empty site (defaults to _old_occ_site)
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S> void pcps::RPF_Descriptor<S>::get_row_col_diff(const std::vector<S> & pm,
                                                                 S * const vec,
                                                                 int new_ind,
                                                                 int old_ind) const {

  // if not supplied, default the indices to those held internally
  new_ind = ( new_ind < 0 ? _new_occ_site : new_ind );
  old_ind = ( old_ind < 0 ? _old_occ_site : old_ind );

  assert( new_ind >= 0 );
  assert( old_ind >= 0 );
  assert( new_ind < 2*_nas );
  assert( old_ind < 2*_nas );

  // get whether the move was alpha or beta
  const int b = old_ind % 2;
  assert( new_ind % 2 == old_ind % 2 );

  // switch to alpha or beta only indexing
  new_ind /= 2;
  old_ind /= 2;

  // compute the row or column difference
  if (b == 0)
    for (int i = 0; i < _nap; i++)
      vec[i] = pm[ new_ind * _nas + _b_sites[i] ] -  pm[ old_ind * _nas + _b_sites[i] ];
  else
    for (int i = 0; i < _nap; i++)
      vec[i] = pm[ _a_sites[i] * _nas + new_ind ] -  pm[ _a_sites[i] * _nas + old_ind ];

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  returns the sign change due to implicit row or column reordering for the move stored
///         by the descriptor
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S> double pcps::RPF_Descriptor<S>::sign_change() const {
  assert( _new_occ_site % 2 == _old_occ_site % 2 );
  const int b = _new_occ_site % 2;
  const int rc_ind = this->get_row_col_index(_old_occ_site);
  const int new_ind = _new_occ_site / 2;
  const int old_ind = _old_occ_site / 2;
  const std::vector<int> & occ_sites = ( b == 0 ? _a_sites : _b_sites );
  const std::vector<int> & ordered_perm = ( b == 0 ? _a_ordered_perm : _b_ordered_perm );
  const std::vector<int> & perm_pos = ( b == 0 ? _a_perm_pos : _b_perm_pos );
  int count = 0;
  for (int i = perm_pos[rc_ind] + 1; i < _nap && occ_sites[ordered_perm[i]] < new_ind; i++, count++) {}
  for (int i = perm_pos[rc_ind] - 1; i >= 0 && occ_sites[ordered_perm[i]] > new_ind; i--, count++) {}
  return 1.0 - 2.0 * ( count % 2 );
}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  returns the sign due to reordering of the current occupation vectors
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S> double pcps::RPF_Descriptor<S>::reordering_sign() const {

  struct SpinSign {
    static double get(const std::vector<int> & unsorted) {
      std::vector<int> v(unsorted);
      double retval = 1.0;
      for (int i = 0; i < v.size()-1; ) {
        if ( v.at(i) > v.at(i+1) ) {
          std::swap(v.at(i), v.at(i+1));
          retval = -retval;
          i = 0;
        } else {
          i++;
        }
      }
      return retval;
    }
  };

  return SpinSign::get(_a_sites) * SpinSign::get(_b_sites);

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief Constructs the pairing function by reading it from an archive.
///
/// \param[in]      userinp  the user's input options
/// \param[in,out]  arch     an archive from which to read the correlator
/// \param[in]      wi_map   map of names to available integer wavefunction data objects
/// \param[in]      wf_map   map of names to available floating point wavefunction data objects
/// \param[in]      co_map   map of names to available coefficient objects
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> pcps::RestPairingFunc<S>::RestPairingFunc(const pcps::Input & userinp,
                                                             pcps::Archive & arch,
                                                             const std::map<std::string, boost::shared_ptr<pcps::WfnData<int> > > & wi_map,
                                                             const std::map<std::string, boost::shared_ptr<pcps::WfnData<S> > > & wf_map,
                                                             const std::map<std::string, boost::shared_ptr<pcps::CoeffObj<S> > > & co_map)
  : CoeffObj<S>(1, 0, 0),
    _nas(userinp.nsites()/2),
    _nap(userinp.nparticles()/2),
    _nau(userinp.nsites()/2 - userinp.nparticles()/2),
    _best_prop_p(-1),
    _best_prop_q(-1),
    _wfn_sampling(true),
    _pqpq_guiding(userinp.pqpq_guiding()),
    _diag_rdm_damping(userinp.diag_rdm_damping()),
    _best_prop_ratio(pcps::unity<S>()),
    _acti_ratio_sum(pcps::zero<S>()),
    _prop_ratio_sum(pcps::zero<S>()),
    _init_value(pcps::unity<S>()),
    _rest_value(pcps::unity<S>()),
    _acti_a_ratios( (userinp.nparticles()/2) * ( userinp.nsites()/2 - userinp.nparticles()/2 + 1 ), pcps::zero<S>()),
    _acti_b_ratios( (userinp.nparticles()/2) * ( userinp.nsites()/2 - userinp.nparticles()/2 + 1 ), pcps::zero<S>()),
    _prop_a_ratios( (userinp.nparticles()/2) * ( userinp.nsites()/2 - userinp.nparticles()/2 + 1 ), pcps::zero<S>()),
    _prop_b_ratios( (userinp.nparticles()/2) * ( userinp.nsites()/2 - userinp.nparticles()/2 + 1 ), pcps::zero<S>()),
    _Zmat(userinp.nparticles() * userinp.nparticles() / 4, pcps::zero<S>()),
    _ZTmat(userinp.nparticles() * userinp.nparticles() / 4, pcps::zero<S>()),
    _Umat(4 * userinp.nparticles() / 2, pcps::zero<S>()),
    _Vmat(4 * userinp.nparticles() / 2, pcps::zero<S>()),
    _VZmat(4 * userinp.nparticles() / 2, pcps::zero<S>()),
    _UZTmat(4 * userinp.nparticles() / 2, pcps::zero<S>()),
    _Qmat(4*4, pcps::zero<S>()),
    _Imat(4*4, pcps::zero<S>()),
    _work(std::max(userinp.nparticles() * userinp.nparticles() / 4, 2*4*4 + 4 * userinp.nparticles() / 2), pcps::zero<S>()),
    _pqpq_tei( userinp.nsites() * userinp.nsites() / 4 ),
    _Umask(4, 0),
    _Vmask(4, 0),
    _Qmask(4*4, 0),
    _iwork(2 * std::max(3, userinp.nparticles()/2), 0),
    _rest_desc(userinp),
    _acti_desc(userinp),
    _prop_desc(userinp)
{

  // check that we have zero net spin
  if (userinp.z_spin() != 0.0)
    throw pcps::Exception("the restricted pairing function is only defined for total Sz = 0");

  // check that we an even particle number
  if (userinp.nparticles() % 2 != 0)
    throw pcps::Exception("the restricted pairing function is only defined for even particle number");

  // get name of the pairing function
  pcps::mpi::read_and_bcast(arch, this->_name, "Failed to read name in RestPairingFunc constructor.");

  // get values
  std::string wf_name;
  pcps::mpi::read_and_bcast(arch, wf_name, "Failed to read name of value data in RestPairingFunc constructor.");
  if (wf_map.count(wf_name) == 0)
    throw pcps::Exception("no floating point data by the name %s") % wf_name;
  this->_wf_deps->at(0) = wf_map.find(wf_name)->second;

  // check that the value data object has the correct number of elements
  if ( this->pairing_matrix().size() != userinp.nsites() * userinp.nsites() / 4 )
    throw pcps::Exception("Wrong number of values in RestPairingFunc constructor.  Found %i values, expected %i.")
          % this->pairing_matrix().size() % ( userinp.nsites() * userinp.nsites() / 4 );

  // initialize the correct value map
  this->init_correct_map();

}

template pcps::RestPairingFunc<double>::RestPairingFunc
(const pcps::Input & userinp,
 pcps::Archive & arch,
 const std::map<std::string, boost::shared_ptr<pcps::WfnData<int> > > & wi_map,
 const std::map<std::string, boost::shared_ptr<pcps::WfnData<double> > > & wf_map,
 const std::map<std::string, boost::shared_ptr<pcps::CoeffObj<double> > > & co_map);

template pcps::RestPairingFunc<std::complex<double> >::RestPairingFunc
(const pcps::Input & userinp,
 pcps::Archive & arch,
 const std::map<std::string, boost::shared_ptr<pcps::WfnData<int> > > & wi_map,
 const std::map<std::string, boost::shared_ptr<pcps::WfnData<std::complex<double> > > > & wf_map,
 const std::map<std::string, boost::shared_ptr<pcps::CoeffObj<std::complex<double> > > > & co_map);

// ///////////////////////////////////////////////////////////////////////////////////////////////////
// /// \brief   computes the products between the inverse of the occpied pairing matrix and the
// ///          various noise vectors
// ///
// /// \param[in]      config  the lattice configuration
// ///
// ///////////////////////////////////////////////////////////////////////////////////////////////////
// template<class S> void pcps::RestPairingFunc<S>::compute_noise_products(const pcps::Configuration<S> & config) {
// 
//   const std::vector<S> & u1 = config.noise_vecs_u().at(0);
//   const std::vector<S> & u2 = config.noise_vecs_u().at(1);
//   const std::vector<S> & v1 = config.noise_vecs_v().at(0);
//   const std::vector<S> & v2 = config.noise_vecs_v().at(1);
// 
//   pcps::xgemm('T', 'N', _nap, 1, _nap, pcps::unity<S>(), &_inv_mat.at(0), _nap, &u1.at(0), _nap, pcps::zero<S>(), &_Iu1.at(0), _nap);
//   pcps::xgemm('T', 'N', _nap, 1, _nap, pcps::unity<S>(), &_inv_mat.at(0), _nap, &u2.at(0), _nap, pcps::zero<S>(), &_Iu2.at(0), _nap);
// 
//   pcps::xgemm('N', 'N', _nap, 1, _nap, pcps::unity<S>(), &_inv_mat.at(0), _nap, &v1.at(0), _nap, pcps::zero<S>(), &_v1I.at(0), _nap);
//   pcps::xgemm('N', 'N', _nap, 1, _nap, pcps::unity<S>(), &_inv_mat.at(0), _nap, &v2.at(0), _nap, pcps::zero<S>(), &_v2I.at(0), _nap);
// 
//   //pcps::xgemm('N', 'T', 1, _nap, _nap, pcps::unity<S>(), &v1.at(0), 1, &_inv_mat.at(0), _nap, pcps::zero<S>(), &_v1I.at(0), 1);
//   //pcps::xgemm('N', 'T', 1, _nap, _nap, pcps::unity<S>(), &v2.at(0), 1, &_inv_mat.at(0), _nap, pcps::zero<S>(), &_v2I.at(0), 1);
// 
// }

// ///////////////////////////////////////////////////////////////////////////////////////////////////
// /// \brief prints out the current occupancy and permutation vectors
// ///
// ///////////////////////////////////////////////////////////////////////////////////////////////////
// template<class S> void pcps::RestPairingFunc<S>::print_occupancy_info() {
// 
//   std::cout << boost::format("row sites     (alpha):");
//   for (int i = 0; i < _nap; i++)
//     std::cout << boost::format(" %3i") % _a_sites.at(i);
//   std::cout << std::endl;
// 
//   std::cout << boost::format("row ordering  (alpha):");
//   for (int i = 0; i < _nap; i++)
//     std::cout << boost::format(" %3i") % _a_ordered_perm.at(i);
//   std::cout << std::endl;
// 
//   std::cout << boost::format("row perm_pos  (alpha):");
//   for (int i = 0; i < _nap; i++)
//     std::cout << boost::format(" %3i") % _a_perm_pos.at(i);
//   std::cout << std::endl;
// 
//   std::cout << std::endl;
// 
//   std::cout << boost::format("col sites     (beta) :");
//   for (int i = 0; i < _nap; i++)
//     std::cout << boost::format(" %3i") % _b_sites.at(i);
//   std::cout << std::endl;
// 
//   std::cout << boost::format("col ordering  (beta) :");
//   for (int i = 0; i < _nap; i++)
//     std::cout << boost::format(" %3i") % _b_ordered_perm.at(i);
//   std::cout << std::endl;
// 
//   std::cout << boost::format("col perm_pos  (beta) :");
//   for (int i = 0; i < _nap; i++)
//     std::cout << boost::format(" %3i") % _b_perm_pos.at(i);
//   std::cout << std::endl;
// 
//   std::cout << std::endl;
// 
// }

// ///////////////////////////////////////////////////////////////////////////////////////////////////
// /// \brief Initializes the object for the given lattice configuration and checks if noise is
// ///        correctly removed using the matrix determinant lemma
// ///
// /// \param[in]      config  the lattice configuration we want the object initialized for
// /// \param[in]      check   whether to check noise removal
// ///
// ///////////////////////////////////////////////////////////////////////////////////////////////////
// template<class S> void pcps::RestPairingFunc<S>::check_initialize(const pcps::Configuration<S> & config, const bool check) {
// 
//   // populate lists of occupied alpha and beta sites
//   this->populate_site_lists(config);
// 
//   // build the occupied pairing matrix whose inverse and determinant we need
//   this->build_occupied_matrix(config, _inv_mat);
// 
//   // compute the inverse matrix and determinant using an LU decomposition
//   S det_val;
//   pcps::matrix_inverse_lu(_nap, det_val, &_inv_mat.at(0), &_work_mat.at(0), &_iwork.at(0));
// 
//   // set the current value
//   this->_curr_val = pcps::MantExp<S>(det_val);
// 
//   // compute the products of the inverse times the noise vectors
//   this->compute_noise_products(config);
// 
//   //this->print_occupancy_info();
// 
//   // test that noise removal is correct
//   if (check) {
//     std::cout << "checking noise removal after initialization:" << std::endl << std::endl;
//     this->test_noise_removal(config);
//   }
//   
// }

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief Get the [pq,pq] two electron integrals from the supplied Hamiltonian
///
/// \param[in]      ham     the Hamiltonian
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S> void pcps::RestPairingFunc<S>::get_pqpq_ints(const pcps::Hamiltonian<S> & ham) {

  ham.get_pqpq_ints(_pqpq_tei);

}

template void pcps::RestPairingFunc<double>::get_pqpq_ints(const pcps::Hamiltonian<double> & ham);
template void pcps::RestPairingFunc<std::complex<double> >::get_pqpq_ints(const pcps::Hamiltonian<std::complex<double> > & ham);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief Initializes the object for the given lattice configuration.
///
/// \param[in]      config  the lattice configuration we want the object initialized for
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S> void pcps::RestPairingFunc<S>::initialize(const pcps::Configuration<S> & config) {

  // designate the current configuration as the resting, active, and excited configuration
  _rest_desc.prepare_occ_info(config);
  _acti_desc.prepare_occ_info(config);
  _prop_desc.prepare_occ_info(config);

  // build the occupied pairing matrix for the resting configuration
  _rest_desc.build_occupied_matrix(this->pairing_matrix(), _Zmat);
  _rest_desc.build_occupied_matrix(this->pairing_matrix(), _ZTmat);

  //for (int i = 0; i < _nap; i++) {
  //  for (int j = 0; j < _nap; j++)
  //    std::cout << boost::format("  %20.12f") % _ZTmat.at(i*_nap+j);
  //  std::cout << std::endl;
  //}
  //std::cout << std::endl;

  // set the current value, the initialization value, and the resting value
  this->_curr_val = pcps::matrix_determinant_lu(_nap, &_ZTmat.at(0), &_iwork.at(0));
  _init_value = this->_curr_val;
  _rest_value = this->_curr_val;

  //std::cout << boost::format("  %20.12f") % this->_curr_val.get() << std::endl;
  //std::cout << std::endl;

  //for (int i = 0; i < _nap; i++) {
  //  for (int j = 0; j < _nap; j++)
  //    std::cout << boost::format("  %20.12f") % _Zmat.at(i*_nap+j);
  //  std::cout << std::endl;
  //}
  //std::cout << std::endl;

  // compute the inverse matrix
  S det_val;
  pcps::matrix_inverse_lu(_nap, det_val, &_Zmat.at(0), &_ZTmat.at(0), &_iwork.at(0));

  // compute the transpose of the inverse matrix
  for (int i = 0; i < _nap; i++)
  for (int j = 0; j < _nap; j++)
    _ZTmat[_nap*i+j] = _Zmat[_nap*j+i];

  // use the get_value function to compute the active configuration's ratios and ratio sum
  {
    pcps::Configuration<S> temp_config(config);
    pcps::ConfigDiff<S> diff;
    diff.newly_full().assign(1, _rest_desc.new_occ_site());
    diff.newly_empty().assign(1, _rest_desc.old_occ_site());
    bool option_holder = false;
    std::swap(_wfn_sampling, option_holder);
    this->get_value(temp_config, config, diff);
    std::swap(_wfn_sampling, option_holder);
    _acti_a_ratios  = _prop_a_ratios;
    _acti_b_ratios  = _prop_b_ratios;
    _acti_ratio_sum = _prop_ratio_sum;
    _acti_desc = _rest_desc;
    _prop_desc = _rest_desc;
  }

}

template void pcps::RestPairingFunc<double>::initialize(const pcps::Configuration<double> &);
template void pcps::RestPairingFunc<std::complex<double> >::initialize(const pcps::Configuration<std::complex<double> > &);

// ///////////////////////////////////////////////////////////////////////////////////////////////////
// /// \brief builds the little 3x3 matrix involved in the matrix determinant lemma and the
// ///        Woodbury matrix identity.
// ///
// /// \param[in]      diff        an object specifying how the configurations differ
// ///
// ///////////////////////////////////////////////////////////////////////////////////////////////////
// template<class S> void pcps::RestPairingFunc<S>::build_little_mat(const pcps::ConfigDiff<S> & diff) {
// 
//   // get whether the move was alpha or beta
//   const int b = diff.newly_empty().at(0) % 2;
// 
//   // get the row or column index of the particle that was destroyed
//   const int rc_ind = this->get_row_col_index(diff.newly_empty().at(0));
// 
//   // get the row or column difference
//   this->get_row_col_diff(diff, &_work_mat.at(0));
// 
//   // set the little matrix equal to the identity
//   std::fill(_little_mat.begin(), _little_mat.end(), pcps::zero<S>());
//   _little_mat[0*3+0] = pcps::unity<S>();
//   _little_mat[1*3+1] = pcps::unity<S>();
//   _little_mat[2*3+2] = pcps::unity<S>();
// 
//   // alpha case
//   if (b == 0) {
// 
//     if ( diff.noise_change_codes().at(0) == 'u' || diff.noise_change_codes().at(1) == 'u' )
//       throw pcps::Exception("pcps::RestPairingFunc<S>::build_little_mat requires noise changes of type 'v' for an alpha move");
// 
//     _little_mat[0*3+0] += pcps::xdot(_nap, &_inv_mat.at(rc_ind), _nap, &_work_mat.at(0), 1);
// 
//     _little_mat[0*3+1] += pcps::xdot(_nap, &_Iu1.at(0), 1, &_work_mat.at(0), 1);
//     _little_mat[0*3+2] += pcps::xdot(_nap, &_Iu2.at(0), 1, &_work_mat.at(0), 1);
// 
//     if ( diff.noise_change_codes().at(0) == 'v' ) {
//       _little_mat[1*3+0] += pcps::xdot(_nap, &_inv_mat.at(rc_ind), _nap, &diff.noise_changes().at(0).at(0), 1);
//       _little_mat[1*3+1] += pcps::xdot(_nap, &_Iu1.at(0), 1, &diff.noise_changes().at(0).at(0), 1);
//       _little_mat[1*3+2] += pcps::xdot(_nap, &_Iu2.at(0), 1, &diff.noise_changes().at(0).at(0), 1);
//     }
// 
//     if ( diff.noise_change_codes().at(1) == 'v' ) {
//       _little_mat[2*3+0] += pcps::xdot(_nap, &_inv_mat.at(rc_ind), _nap, &diff.noise_changes().at(1).at(0), 1);
//       _little_mat[2*3+1] += pcps::xdot(_nap, &_Iu1.at(0), 1, &diff.noise_changes().at(1).at(0), 1);
//       _little_mat[2*3+2] += pcps::xdot(_nap, &_Iu2.at(0), 1, &diff.noise_changes().at(1).at(0), 1);
//     }
// 
//   // beta case
//   } else {
// 
//     if ( diff.noise_change_codes().at(0) == 'v' || diff.noise_change_codes().at(1) == 'v' )
//       throw pcps::Exception("pcps::RestPairingFunc<S>::build_little_mat requires noise changes of type 'u' for a beta move");
// 
//     _little_mat[0*3+0] += pcps::xdot(_nap, &_inv_mat.at(rc_ind*_nap), 1, &_work_mat.at(0), 1);
// 
//     _little_mat[1*3+0] += pcps::xdot(_nap, &_v1I.at(0), 1, &_work_mat.at(0), 1);
//     _little_mat[2*3+0] += pcps::xdot(_nap, &_v2I.at(0), 1, &_work_mat.at(0), 1);
// 
//     if ( diff.noise_change_codes().at(0) == 'u' ) {
//       _little_mat[0*3+1] += pcps::xdot(_nap, &_inv_mat.at(rc_ind*_nap), 1, &diff.noise_changes().at(0).at(0), 1);
//       _little_mat[1*3+1] += pcps::xdot(_nap, &_v1I.at(0), 1, &diff.noise_changes().at(0).at(0), 1);
//       _little_mat[2*3+1] += pcps::xdot(_nap, &_v2I.at(0), 1, &diff.noise_changes().at(0).at(0), 1);
//     }
// 
//     if ( diff.noise_change_codes().at(1) == 'u' ) {
//       _little_mat[0*3+2] += pcps::xdot(_nap, &_inv_mat.at(rc_ind*_nap), 1, &diff.noise_changes().at(1).at(0), 1);
//       _little_mat[1*3+2] += pcps::xdot(_nap, &_v1I.at(0), 1, &diff.noise_changes().at(1).at(0), 1);
//       _little_mat[2*3+2] += pcps::xdot(_nap, &_v2I.at(0), 1, &diff.noise_changes().at(1).at(0), 1);
//     }
// 
//   }
// 
// }

// ///////////////////////////////////////////////////////////////////////////////////////////////////
// /// \brief Given the corresponding little matrix, returns the value of the coefficient for a
// ///        configuration that differs from the current configuration by a single move.
// ///
// /// \param[in]      new_ind     the index of the newly occupied site
// /// \param[in]      old_ind     the index of the newly empty site
// /// \param[in]      m           the little matrix, which should already have been computed
// ///
// /// \return the coefficient at the new configuration
// ///
// ///////////////////////////////////////////////////////////////////////////////////////////////////
// template<class S> pcps::MantExp<S> pcps::RestPairingFunc<S>::get_single_excite_value(int new_ind,
//                                                                                      int old_ind,
//                                                                                      const std::vector<S> & m) const
// {
// 
//   // get whether the move was alpha or beta
//   const int b = new_ind % 2;
// 
//   // get the row or column index of the particle that was destroyed
//   const int rc_ind = this->get_row_col_index(old_ind);
// 
//   // change from all site indexing to alpha/beta only site indexing
//   new_ind /= 2;
//   old_ind /= 2;
// 
//   // compute the sign associated with the excitation due to possible changes in the permutation
//   // that orders the rows or columns
//   const std::vector<int> & occ_sites = ( b == 0 ? _a_sites : _b_sites );
//   const std::vector<int> & ordered_perm = ( b == 0 ? _a_ordered_perm : _b_ordered_perm );
//   const std::vector<int> & perm_pos = ( b == 0 ? _a_perm_pos : _b_perm_pos );
//   int count = 0;
//   for (int i = perm_pos[rc_ind] + 1; i < _nap && occ_sites[ordered_perm[i]] < new_ind; i++, count++) {}
//   for (int i = perm_pos[rc_ind] - 1; i >= 0 && occ_sites[ordered_perm[i]] > new_ind; i--, count++) {}
//   const double sign = 1.0 - 2.0 * ( count % 2 );
// 
//   //std::cout << boost::format("sign = %6.1f") % sign << std::endl;
// 
//   // return the coefficient value at the excited configuration, which is the current value times the determinant
//   // of the little matrix times the sign due to row/column reordering
//   assert( m.size() == 9 );
//   return this->_curr_val * ( sign * (   m[0*3+0] * ( m[1*3+1] * m[2*3+2] - m[2*3+1] * m[1*3+2] )
//                                       - m[0*3+1] * ( m[1*3+0] * m[2*3+2] - m[2*3+0] * m[1*3+2] )
//                                       + m[0*3+2] * ( m[1*3+0] * m[2*3+1] - m[2*3+0] * m[1*3+1] ) ) );
// 
// }

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief Populates the V or U vector corresponding to a single electron move
///        described by the supplied descriptor.
///
/// \param[in]      vp       which vector to populate (the first, the second, etc.)
/// \param[in]      b        0-alpha, 1-beta
/// \param[in]      d        the object describing the move
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S> void pcps::RestPairingFunc<S>::get_update_vec_V_U(const int vp,
                                                                    const int b,
                                                                    const pcps::RPF_Descriptor<S> & d,
                                                                    const int rc_ind) {

  S * const V = ( b == 0 ? &_Vmat[0] : &_Umat[0] );
  d.get_row_col_diff(this->pairing_matrix(), V + _nap*vp);
  ( b == 0 ? _Vmask[vp] : _Umask[vp] ) = 0;

  S * const U = ( b == 0 ? &_Umat[0] : &_Vmat[0] );
  std::fill(U + vp*_nap, U + (vp+1)*_nap, pcps::zero<S>());
  U[vp*_nap + rc_ind] = pcps::unity<S>();

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief Populates the UZT or VZ vector corresponding to a single electron move in which
///        the destroyed particle corresponds to the given row or column.
///
/// \param[in]      vp       which vector to populate (the first, the second, etc.)
/// \param[in]      b        0-alpha, 1-beta
/// \param[in]      rc_ind   the row or column of the destroyed electron
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S> void pcps::RestPairingFunc<S>::get_update_vec_UZT_VZ(const int vp,
                                                                       const int b,
                                                                       const int rc_ind) {

  const S * const Z = ( b == 0 ? &_Zmat[0] : &_ZTmat[0] );
  S * const UZT = ( b == 0 ? &_UZTmat[0] : &_VZmat[0] );
  pcps::xcopy(_nap, Z + rc_ind, _nap, UZT + _nap*vp, 1);
  ( b == 0 ? _Umask[vp] : _Vmask[vp] ) = 1;

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief Populates the V and UZT or U and VZ vectors corresponding to the single electron move
///        described by the supplied descriptor.
///
/// \param[in]      vp       which vectors to populate (the first set, the second set, etc.)
/// \param[in]      d        the object describing the move
/// \param[in]      rc_ind   the row or column of the destroyed electron (computed on the fly
///                          if a negative value is given)
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S> void pcps::RestPairingFunc<S>::get_update_vecs(const int vp,
                                                                 const pcps::RPF_Descriptor<S> & d,
                                                                 int rc_ind) {

  // get the spin for the move
  const int b = d.new_occ_site() % 2;

  // if not given or given as a negative value, get the row/column index of the destroyed particle for the move
  if (rc_ind < 0)
    rc_ind = d.get_row_col_index(d.old_occ_site());

  // get V or U vector
  this->get_update_vec_V_U(vp, b, d, rc_ind);

  // get UZT or VZ vector
  this->get_update_vec_UZT_VZ(vp, b, rc_ind);

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief returns the value of the coefficient for a configuration that differs from the
///        currently initialized configuration as specified.
///
/// \param[in]      new_config  the lattice configuration whose coefficient we want
/// \param[in]      old_config  the lattice configuration the object is currently initialized for
/// \param[in]      diff        an object specifying how the configurations differ
///
/// \return the coefficient at the new configuration
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S> pcps::MantExp<S> pcps::RestPairingFunc<S>::get_value(const pcps::Configuration<S> & new_config,
                                                                       const pcps::Configuration<S> & old_config,
                                                                       const pcps::ConfigDiff<S> & diff)
{

  pcps::start_timer("get_value");

  // check that exactly one particle moved
  if ( diff.newly_empty().size() != 1 || diff.newly_full().size() != 1 )
    throw pcps::Exception("pcps::RestPairingFunc<S>::get_value currently only supports single particle moves");

  const int new_ind = diff.newly_full().at(0);
  const int old_ind = diff.newly_empty().at(0);

  // check that the move did not change the particle's spin
  if ( diff.newly_empty().at(0) % 2 != diff.newly_full().at(0) % 2 )
    throw pcps::Exception("pcps::RestPairingFunc<S>::get_value requires the two sites to be of the same spin");

  // get a reference to the fock configuration vector for the proposed configuration
  const std::vector<int> & fc = new_config.fock_config();

  // initialize the return value
  pcps::MantExp<S> retval;

  // if we are sampling with the wavefunction as the guiding function
  if ( _wfn_sampling ) {

    // get whether the excitation is alpha or beta
    const int b = new_ind % 2;

    if (true) {

      // get the U and V vectors for the single excitation
      _prop_desc = _rest_desc;
      _prop_desc.new_occ_site() = new_ind;
      _prop_desc.old_occ_site() = old_ind;
      this->get_update_vecs(0, _prop_desc);

      // compute the value at the excited configuration
      _Imat[0] = pcps::unity<S>();
      _Qmask[0] = 0;
      retval = this->_curr_val * pcps::matrix_determinant_lemma(_nap, 1, &_Zmat[0], &_Umat[0], &_Vmat[0], &_Imat[0],
                                                                &_UZTmat[0], &_VZmat[0], &_Qmat[0],
                                                                &_Umask[0], &_Vmask[0], &_Qmask[0], &_work[0], &_iwork[0], true);

    } else {

      std::fill(_Imat.begin(), _Imat.end(), pcps::zero<S>());
      _Imat[0*3+0] = pcps::unity<S>();
      _Imat[1*3+1] = pcps::unity<S>();
      _Imat[2*3+2] = pcps::unity<S>();

      std::fill(_Umat.begin(), _Umat.end(), pcps::zero<S>());
      std::fill(_Vmat.begin(), _Vmat.end(), pcps::zero<S>());
      std::fill(_VZmat.begin(), _VZmat.end(), pcps::zero<S>());
      std::fill(_UZTmat.begin(), _UZTmat.end(), pcps::zero<S>());
      std::fill(_Qmat.begin(), _Qmat.end(), pcps::zero<S>());

      std::fill(_Umask.begin(), _Umask.end(), 0);
      std::fill(_Vmask.begin(), _Vmask.end(), 0);
      std::fill(_Qmask.begin(), _Qmask.end(), 0);

      _prop_desc = _rest_desc;
      _prop_desc.new_occ_site() = new_ind;
      _prop_desc.old_occ_site() = old_ind;
      this->get_update_vecs(2, _prop_desc);

      retval = this->_curr_val * pcps::matrix_determinant_lemma(_nap, 3, &_Zmat[0], &_Umat[0], &_Vmat[0], &_Imat[0],
                                                                &_UZTmat[0], &_VZmat[0], &_Qmat[0],
                                                                &_Umask[0], &_Vmask[0], &_Qmask[0], &_work[0], &_iwork[0], true);

    }

    // get any sign change due to implicit row/column reordering
    retval *= _prop_desc.sign_change();

  // if we are sampling with the local rdm as the guiding function
  } else {

    //pcps::start_timer("prep");

    // prepare the identity matrix for a rank three update
    std::fill(_Imat.begin(), _Imat.end(), pcps::zero<S>());
    _Imat[0*3+0] = pcps::unity<S>();
    _Imat[1*3+1] = pcps::unity<S>();
    _Imat[2*3+2] = pcps::unity<S>();

    // reset the masks
    std::fill(_Umask.begin(), _Umask.end(), 0);
    std::fill(_Vmask.begin(), _Vmask.end(), 0);
    std::fill(_Qmask.begin(), _Qmask.end(), 0);

    // prepare information for the active to proposed move
    _acti_desc.new_occ_site() = new_ind;
    _acti_desc.old_occ_site() = old_ind;

    // get the sign change from implicit reordering in rest->acti and acti->prop moves
    double base_sign = _rest_desc.sign_change();
    base_sign *= _acti_desc.sign_change();

    // get a description of the pairing function at the proposed configuration
    _prop_desc = _acti_desc;
    _prop_desc.update(new_ind, old_ind);

    // get the update vectors for the move from resting to active
    this->get_update_vecs(0, _rest_desc);

    // get the update vectors for the move from active to proposed
    this->get_update_vecs(1, _acti_desc);

    //// print some info
    //if (false) {
    //  std::cout << boost::format("%6s  %16s  %16s  %16s  %16s  %16s  %16s  %16s")
    //               % ""
    //               % "a_sites"
    //               % "b_sites"
    //               % "a_perm"
    //               % "b_perm"
    //               % "a_ppos"
    //               % "b_ppos"
    //               % "move"
    //           << std::endl;
    //  _rest_desc.printout("rest");
    //  _acti_desc.printout("acti");
    //  _prop_desc.printout("prop");
    //  //std::cout << std::endl;
    //}

    // get the spins for the rest->acti and acti->prop moves
    const int rest_b = _rest_desc.new_occ_site() % 2;
    const int acti_b = _acti_desc.new_occ_site() % 2;

    // get the row or column indices for the rest->acti and acti->prop moves
    const int rest_rc_ind = _rest_desc.get_row_col_index(_rest_desc.old_occ_site());
    const int acti_rc_ind = _acti_desc.get_row_col_index(_acti_desc.old_occ_site());

    // initialize the sum of square ratios
    _prop_ratio_sum = pcps::zero<S>();

    // get the ratio of the resting to the initial determinant values
    const S rest_init_ratio = (_rest_value/_init_value).get();

    // initialize the indices for the best next resting function
    _best_prop_p = -1;
    _best_prop_q = -1;
    double biggest_ratio = 0.0;

    //pcps::stop_timer("prep");

    // get the ratios   < n | p^(+) q | Psi(proposed) > / < n | Psi(initialized) >
    for (int b = 0; b < 2; b++) {

      //pcps::start_timer("prep");

      // get the vector to place the ratios in
      std::vector<S> & prop_ratios = ( b == 0 ? _prop_a_ratios : _prop_b_ratios );

      // zero the ratios
      std::fill(prop_ratios.begin(), prop_ratios.end(), pcps::zero<S>());

      // initialize the ratio storage row counter
      int p_count = 0;

      //pcps::stop_timer("prep");

      //std::cout << "b = " << b << std::endl;

      // loop over the p operator
      for (int p = 0; p < _nas; p++) {

        // skip excitations that destroy the bra
        if ( fc[2*p+b] == 0 ) continue;

        //pcps::start_timer("prep");

        // get the row or column index of the particle that was destroyed
        const int rc_ind = _prop_desc.get_row_col_index(2*p+b);

        // get the UZT or VZ vector for the proposed to excited move
        this->get_update_vec_UZT_VZ(2, b, rc_ind);

        // reset the appropriate Qmask elements
        _Qmask[ 0 * (3-2*b) + 2 * (1+2*b) ] = 0;
        _Qmask[ 1 * (3-2*b) + 2 * (1+2*b) ] = 0;
        _Qmask[ 2 * (3-2*b) + 2 * (1+2*b) ] = 0;

        // initialize the ratio storage column counter
        int q_count = 0;

        //pcps::stop_timer("prep");

        // loop over the q operator
        for (int q = 0; q < _nas; q++) {

          // skip excitations that destroy the bra
          if ( q != p && fc[2*q+b] == 1 ) continue;

          //pcps::start_timer("prep");

          // prepare information for the proposed to excited move
          _prop_desc.new_occ_site() = 2*q+b;
          _prop_desc.old_occ_site() = 2*p+b;

          // get the V or U vector
          this->get_update_vec_V_U(2, b, _prop_desc, rc_ind);

          //this->get_update_vecs(2, _prop_desc);

          // reset the appropriate Qmask elements
          _Qmask[ 2 * (3-2*b) + 0 * (1+2*b) ] = 0;
          _Qmask[ 2 * (3-2*b) + 1 * (1+2*b) ] = 0;
          _Qmask[ 2 * (3-2*b) + 2 * (1+2*b) ] = 0;

          // compute the Q elements in which both the U and V vectors are elementary
          if ( rest_b == 0 && acti_b == 1) { _Qmat[1*3+0] = _Zmat[ acti_rc_ind * _nap + rest_rc_ind ]; _Qmask[1*3+0] = 1; }
          if ( rest_b == 1 && acti_b == 0) { _Qmat[0*3+1] = _Zmat[ rest_rc_ind * _nap + acti_rc_ind ]; _Qmask[0*3+1] = 1; }
          if ( rest_b == 0 &&      b == 1) { _Qmat[2*3+0] = _Zmat[      rc_ind * _nap + rest_rc_ind ]; _Qmask[2*3+0] = 1; }
          if ( rest_b == 1 &&      b == 0) { _Qmat[0*3+2] = _Zmat[ rest_rc_ind * _nap +      rc_ind ]; _Qmask[0*3+2] = 1; }
          if ( acti_b == 0 &&      b == 1) { _Qmat[2*3+1] = _Zmat[      rc_ind * _nap + acti_rc_ind ]; _Qmask[2*3+1] = 1; }
          if ( acti_b == 1 &&      b == 0) { _Qmat[1*3+2] = _Zmat[ acti_rc_ind * _nap +      rc_ind ]; _Qmask[1*3+2] = 1; }

          //std::fill(_Qmask.begin(), _Qmask.end(), 0);
          //std::fill(_Umask.begin(), _Umask.end(), 0);
          //std::fill(_Vmask.begin(), _Vmask.end(), 0);

          //pcps::stop_timer("prep");
          //pcps::start_timer("det");

          // compute the ratio
          S ratio = pcps::matrix_determinant_lemma(_nap, 3, &_Zmat[0],
                                                   &_Umat[0], &_Vmat[0], &_Imat[0],
                                                   &_UZTmat[0], &_VZmat[0], &_Qmat[0],
                                                   &_Umask[0], &_Vmask[0], &_Qmask[0],
                                                   &_work[0], &_iwork[0], true);

          //pcps::stop_timer("det");
          //pcps::start_timer("sign");

          // apply the signs associated with implicit row/column reorderings
          ratio *= base_sign * _prop_desc.sign_change();

          //// check the sign
          //if (true) {
          //  pcps::RPF_Descriptor<S> exci_desc(_prop_desc);
          //  exci_desc.update(_prop_desc.new_occ_site(), _prop_desc.old_occ_site());
          //  const double sign_prod = base_sign * _prop_desc.sign_change();
          //  const double abs_sign = _rest_desc.reordering_sign() * exci_desc.reordering_sign();
          //  assert( abs_sign == sign_prod );
          //}

          // apply the sign associated with the excitation operator
          {
            const int min_pos = 2 * std::min(p,q) + b + 2;
            const int max_pos = 2 * std::max(p,q) + b;
            int count = 0;
            for (int i = min_pos; i < max_pos; i+=2)
              count += fc[i];
            if ( count % 2 != 0 )
              ratio = -ratio;
            //assert( double(count % 2) == _acti_desc.sign_change() * _prop_desc.sign_change() );
            //const int min_pos = 2 * std::min(p,q) + b + 1;
            //const int max_pos = 2 * std::max(p,q) + b;
            //if ( std::accumulate(config.fock_config().begin()+min_pos, config.fock_config().begin()+max_pos, 0) % 2 == 1)
            //  ratio = -ratio;
          }

          //pcps::stop_timer("sign");

          // save the best excitation for possible use as the next resting configuration
          if ( pcps::square_norm(ratio) > biggest_ratio ) {
            _best_prop_ratio = ratio;
            _best_prop_p = 2*p+b;
            _best_prop_q = 2*q+b;
            biggest_ratio = pcps::square_norm(ratio);
          }

          // We currently have the ratio < n | p^(+) q | Psi(proposed) > / < n | Psi(resting) >.
          // Now convert to the ratio < n | p^(+) q | Psi(proposed) > / < n | Psi(initialized) >.
          ratio *= rest_init_ratio;

          // add the square ratio to the sum of the squares
          const S pqpq_factor = ( _pqpq_guiding ? 0.5 * _pqpq_tei[_nas*p+q] : pcps::unity<S>() );
          const double damping_factor = 1.0 - ( p == q ? _diag_rdm_damping : 0.0 );
          _prop_ratio_sum += pcps::square_norm(ratio) * pqpq_factor * damping_factor;

          // store the ratio
          //prop_ratios[ (_nau+1) * p_count + q_count ] = ratio;
          prop_ratios[ new_config.occ_unocc_cmpd(2*p+b, 2*q+b) ] = ratio;

          //std::cout << new_config.occ_unocc_cmpd(2*p+b, 2*q+b) << std::endl;

          // increment the ratio storage column counter
          q_count++;

        }

        // increment the ratio storage row counter
        p_count++;

      }

      //std::cout << std::endl;

//      // We currently have the ratios < n | p^(+) q | Psi(proposed) > / < n | Psi(resting) >.
//      // Now convert to the ratios < n | p^(+) q | Psi(proposed) > / < n | Psi(initialized) >.
//      //pcps::xscal(prop_ratios.size(), (_rest_value).get(), &prop_ratios[0], 1);
//      pcps::xscal(prop_ratios.size(), (_rest_value/_init_value).get(), &prop_ratios[0], 1);
//
//      // add the sum of the squares of the ratios to the total
//      _prop_ratio_sum += pcps::xdotc(prop_ratios.size(), &prop_ratios[0], 1, &prop_ratios[0], 1);
//
//      // get a reference to the vector of occupied sites
//      const std::vector<int> & occ_to_site   = ( b == 0 ? new_config.a_occ_to_site()   : new_config.b_occ_to_site()   );
//
//      // remove diagonal terms from the sum of squares according to their damping factor
//      for (int i = 0; i < occ_to_site.size(); i++) {
//        const int p = occ_to_site[i];
//        const int cpq = new_config.occ_unocc_cmpd(2*p+b, 2*p+b);
//        _prop_ratio_sum -= _diag_rdm_damping * prop_ratios[cpq] * prop_ratios[cpq];
//      }

    }

    // check the validity of the damping
    if ( _diag_rdm_damping < 0.0 || _diag_rdm_damping > 1.0 )
      throw pcps::Exception("diagonal rdm damping must be between zero (no damping) and one (completely damped)");

    // divide the sum by the normalization constant (the rdm trace)
//    _prop_ratio_sum /= 2.0 * double(_nap) * ( double(_nau) + 1.0 - _diag_rdm_damping );
    if (!_pqpq_guiding)
      _prop_ratio_sum /= 2.0 * double(_nap) * ( double(_nau) + 1.0 - _diag_rdm_damping );

    // set the return value equal to the sum
    retval = pcps::MantExp<S>(_prop_ratio_sum);

  }

  pcps::stop_timer("get_value");

  return retval;

}

template pcps::MantExp<double> pcps::RestPairingFunc<double>::get_value(const pcps::Configuration<double> &,
                                                                        const pcps::Configuration<double> &,
                                                                        const pcps::ConfigDiff<double> &);

template pcps::MantExp<std::complex<double> >
  pcps::RestPairingFunc<std::complex<double> >::get_value(const pcps::Configuration<std::complex<double> > &,
                                                          const pcps::Configuration<std::complex<double> > &,
                                                          const pcps::ConfigDiff<std::complex<double> > &);

// ///////////////////////////////////////////////////////////////////////////////////////////////////
// /// \brief  For the specified spin, gets the ratios for the single excitations of that spin.
// ///
// /// \param[in]      config      the lattice configuration whose coefficient we want
// /// \param[in]      b           specifies the spin, 0-alpha, 1-beta
// /// \param[in]      a_ratios    a vector in which to place the ratios
// ///
// ///////////////////////////////////////////////////////////////////////////////////////////////////
// template<class S> void pcps::RestPairingFunc<S>::get_single_excitation_a_ratios(const pcps::Configuration<S> & config,
//                                                                                 const int b,
//                                                                                 std::vector<S> & a_ratios)
// {
// 
//   // prepare configuration object to represent the excited configuration
//   _work_config = config;
//   if ( b == 0 ) {
//     std::fill(_work_config.noise_vecs_v().at(0).begin(), _work_config.noise_vecs_v().at(0).end(), pcps::zero<S>());
//     std::fill(_work_config.noise_vecs_v().at(1).begin(), _work_config.noise_vecs_v().at(1).end(), pcps::zero<S>());
//   } else {
//     std::fill(_work_config.noise_vecs_u().at(0).begin(), _work_config.noise_vecs_u().at(0).end(), pcps::zero<S>());
//     std::fill(_work_config.noise_vecs_u().at(1).begin(), _work_config.noise_vecs_u().at(1).end(), pcps::zero<S>());
//   }
// 
//   // prepare the object that will represent the difference between the configurations
//   _work_diff.newly_empty().resize(1);
//   _work_diff.newly_full().resize(1);
//   if ( b == 0 ) {
//     std::fill(_work_diff.noise_change_codes().begin(), _work_diff.noise_change_codes().end(), 'v');
//     for (int i = 0; i < _work_diff.noise_changes().at(0).size(); i++)
//       _work_diff.noise_changes().at(0).at(i) = -config.noise_vecs_v().at(0).at(i);
//     for (int i = 0; i < _work_diff.noise_changes().at(1).size(); i++)
//       _work_diff.noise_changes().at(1).at(i) = -config.noise_vecs_v().at(1).at(i);
//   } else {
//     std::fill(_work_diff.noise_change_codes().begin(), _work_diff.noise_change_codes().end(), 'u');
//     for (int i = 0; i < _work_diff.noise_changes().at(0).size(); i++)
//       _work_diff.noise_changes().at(0).at(i) = -config.noise_vecs_u().at(0).at(i);
//     for (int i = 0; i < _work_diff.noise_changes().at(1).size(); i++)
//       _work_diff.noise_changes().at(1).at(i) = -config.noise_vecs_u().at(1).at(i);
//   }
// 
//   // check that the ratio matrix is the right size
//   if ( a_ratios.size() != _nas * _nas )
//     throw pcps::Exception("ratio matrix is the wrong size in pcps::RestPairingFunc<S>::get_single_excitation_a_ratios");
// 
//   // zero the ratio matrix
//   std::fill(a_ratios.begin(), a_ratios.end(), pcps::zero<S>());
// 
//   // fill in the parts of U and V corresponding to noise removal
//   std::vector<S> U_vec(_nap*3, pcps::zero<S>());
//   std::vector<S> V_vec(_nap*3, pcps::zero<S>());
//   pcps::xaxpy(_nap, -pcps::unity<S>(), &config.noise_vecs_u().at(0).at(0), 1, &U_vec.at(_nap*1), 1);
//   pcps::xaxpy(_nap, -pcps::unity<S>(), &config.noise_vecs_u().at(1).at(0), 1, &U_vec.at(_nap*2), 1);
//   pcps::xaxpy(_nap,  pcps::unity<S>(), &config.noise_vecs_v().at(0).at(0), 1, &V_vec.at(_nap*1), 1);
//   pcps::xaxpy(_nap,  pcps::unity<S>(), &config.noise_vecs_v().at(1).at(0), 1, &V_vec.at(_nap*2), 1);
// 
//   // prepare mask arrays
//   std::vector<int> Umask_vec(3, 0);
//   std::vector<int> Vmask_vec(3, 0);
//   std::vector<int> Qmask_vec(9, 0);
//   int * const Umask = ( b == 0 ? &Umask_vec.at(0) : &Vmask_vec.at(0) );
//   int * const Vmask = ( b == 0 ? &Vmask_vec.at(0) : &Umask_vec.at(0) );
//   int * const Qmask = &Qmask_vec.at(0);
// 
//   // get pointers to U and V matrices
//   S * const U = ( b == 0 ? &U_vec.at(0) : &V_vec.at(0) );
//   S * const V = ( b == 0 ? &V_vec.at(0) : &U_vec.at(0) );
// 
//   // get pointer to the inverse matrix
//   std::vector<S> Z_vec(_nap*_nap, pcps::zero<S>());
//   if (b == 0)
//     for (int i = 0; i < _nap; i++)
//     for (int j = 0; j < _nap; j++)
//       Z_vec[i*_nap+j] = _inv_mat[i*_nap+j];
//   else
//     for (int i = 0; i < _nap; i++)
//     for (int j = 0; j < _nap; j++)
//       Z_vec[i*_nap+j] = _inv_mat[j*_nap+i];
//   S * const Z = &Z_vec.at(0);
// 
//   // get pointers to VZ and UZT
//   std::vector<S> UZT_vec(_nap*3, pcps::zero<S>());
//   std::vector<S> VZ_vec(_nap*3, pcps::zero<S>());
//   S * const UZT = &UZT_vec.at(0);
//   S * const VZ = &VZ_vec.at(0);
// 
//   // get the identity matrix
//   std::vector<S> I_vec(9, pcps::zero<S>());
//   I_vec[0*3+0] = pcps::unity<S>();
//   I_vec[1*3+1] = pcps::unity<S>();
//   I_vec[2*3+2] = pcps::unity<S>();
//   S * const I = &I_vec.at(0);
// 
//   // get space for the Q matrix
//   std::vector<S> Q_vec(9, pcps::zero<S>());
//   S * const Q = &Q_vec.at(0);
// 
//   // get work arrays
//   std::vector<S> work_vec(2*3*3 + _nap*3, pcps::zero<S>());
//   S * const work = &work_vec.at(0);
//   std::vector<int> iwork_vec(2*3, 0);
//   int * const iwork = &iwork_vec.at(0);
// 
//   // get UZT rows corresponding to noise removal
//   pcps::xgemm('N', 'N', 1, _nap, _nap, pcps::unity<S>(), U+_nap*1, 1, Z, _nap, pcps::zero<S>(), UZT+_nap*1, 1);
//   pcps::xgemm('N', 'N', 1, _nap, _nap, pcps::unity<S>(), U+_nap*2, 1, Z, _nap, pcps::zero<S>(), UZT+_nap*2, 1);
//   Umask[1] = 1;
//   Umask[2] = 1;
// 
//   // get the ratios   < n | p^(+) q | Psi(no noise) > / < n | Psi(with noise) >
//   for (int p = 0; p < _nas; p++) {
// 
//     // skip excitations that destroy the bra
//     if (config.fock_config().at(2*p+b) == 0) continue;
// 
//     // get the row or column index of the particle that was destroyed
//     const int rc_ind = this->get_row_col_index(2*p+b);
// 
//     // get UZT row corresponding to row replacement
//     pcps::xcopy(_nap, Z + rc_ind, _nap, UZT + _nap*0, 1);
//     Umask[0] = 1;
// 
//     // reset Qmask for new destroyed particle
//     Qmask[0*3+0] = 0;
//     Qmask[1*3+0] = 0;
//     Qmask[2*3+0] = 0;
// 
//     for (int q = 0; q < _nas; q++) {
// 
//       // diagonal term
//       if (p == q) {
// 
//         a_ratios[_nas*p+q] = ( this->noiseless_value(config) / this->_curr_val ).get();
// 
//       // off-diagonal terms
//       } else if (config.fock_config().at(2*q+b) == 0) {
// 
//         // update the new configuration
//         _work_config.fock_config().at(2*p+b) = 0;
//         _work_config.fock_config().at(2*q+b) = 1;
// 
//         // update the configuration difference
//         _work_diff.newly_empty().at(0) = 2*p+b;
//         _work_diff.newly_full().at(0) = 2*q+b;
// 
//         // get the row/col difference
//         this->get_row_col_diff(_work_diff, V);
// 
//         // reset Vmask and Qmask for new row/col difference
//         Vmask[0] = 0;
//         Qmask[0*3+0] = 0;
//         Qmask[0*3+1] = 0;
//         Qmask[0*3+2] = 0;
// 
//         //std::fill(Vmask, Vmask+3, 0);
//         //std::fill(Qmask, Qmask+9, 0);
//         //Qmask[0*3+0] = 0;
//         //Qmask[1*3+0] = 0;
//         //Qmask[2*3+0] = 0;
//         //Qmask[0*3+1] = 0;
//         //Qmask[0*3+2] = 0;
// 
//         // get the ratio up to the sign
//         S ratio = pcps::matrix_determinant_lemma(_nap, 3, Z, U, V, I, UZT, VZ, Q, Umask, Vmask, Qmask, work, iwork);
// 
//         //this->build_little_mat(_work_diff);
//         //_little_mat[0*3+0] -= pcps::unity<S>();
//         //_little_mat[1*3+1] -= pcps::unity<S>();
//         //_little_mat[2*3+2] -= pcps::unity<S>();
//         //std::fill(Qmask, Qmask+9, 1);
//         //S ratio = pcps::matrix_determinant_lemma(_nap, 3, Z, U, V, I, UZT, VZ, &_little_mat.at(0), Umask, Vmask, Qmask, work, iwork);
// 
//         // compute the sign associated with the excitation due to possible changes in the permutation
//         // that orders the rows or columns
//         {
//           const int new_ind = q;
//           const int old_ind = p;
//           const std::vector<int> & occ_sites = ( b == 0 ? _a_sites : _b_sites );
//           const std::vector<int> & ordered_perm = ( b == 0 ? _a_ordered_perm : _b_ordered_perm );
//           const std::vector<int> & perm_pos = ( b == 0 ? _a_perm_pos : _b_perm_pos );
//           int count = 0;
//           for (int i = perm_pos[rc_ind] + 1; i < _nap && occ_sites[ordered_perm[i]] < new_ind; i++, count++) {}
//           for (int i = perm_pos[rc_ind] - 1; i >= 0 && occ_sites[ordered_perm[i]] > new_ind; i--, count++) {}
//           const double sign = 1.0 - 2.0 * ( count % 2 );
//           ratio *= sign;
//         }
// 
//         //// get the ratio
//         //S ratio = ( this->get_value(_work_config, config, _work_diff) / this->_curr_val ).get();
// 
//         // reset the work configuration
//         _work_config.fock_config().at(2*p+b) = 1;
//         _work_config.fock_config().at(2*q+b) = 0;
// 
//         // compute the sign
//         {
//           const int min_pos = 2 * std::min(p,q) + b + 2;
//           const int max_pos = 2 * std::max(p,q) + b;
//           int count = 0;
//           for (int i = min_pos; i < max_pos; i+=2)
//             count += config.fock_config()[i];
//           if ( count % 2 != 0 )
//             ratio = -ratio;
//           //const int min_pos = 2 * std::min(p,q) + b + 1;
//           //const int max_pos = 2 * std::max(p,q) + b;
//           //if ( std::accumulate(config.fock_config().begin()+min_pos, config.fock_config().begin()+max_pos, 0) % 2 == 1)
//           //  ratio = -ratio;
//         }
// 
//         // store the ratio
//         a_ratios[_nas*p+q] = ratio;
// 
//       }
// 
//     }
// 
//   }
// 
// }

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief Computes the ratios of the noiseless coefficients at each singly excited configuration
///        to the noisy coefficient at the current configuration
///
/// \param[in]      config      the lattice configuration the object is currently initialized for
/// \param[out]     a_ratios    on exit, the matrix of excitation ratios for alpha excitations
/// \param[out]     b_ratios    on exit, the matrix of excitation ratios for beta excitations
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S> void pcps::RestPairingFunc<S>::single_excitation_ratios(const pcps::Configuration<S> & config,
                                                                          std::vector<S> & a_ratios,
                                                                          std::vector<S> & b_ratios)
{

  //this->get_single_excitation_a_ratios(config, 0, a_ratios);
  //this->get_single_excitation_a_ratios(config, 1, b_ratios);

  // if we are sampling with the wavefunction as the guiding function
  if ( _wfn_sampling ) {

    // use the value computation function for sampling via the local rdm to compute the ratios
    pcps::Configuration<S> temp_config(config);
    pcps::ConfigDiff<S> diff;
    diff.newly_full().assign(1, _rest_desc.new_occ_site());
    diff.newly_empty().assign(1, _rest_desc.old_occ_site());
    _wfn_sampling = false;
    this->get_value(config, temp_config, diff);
    _wfn_sampling = true;
    a_ratios = _prop_a_ratios;
    b_ratios = _prop_b_ratios;
    _acti_desc = _rest_desc;
    _prop_desc = _rest_desc;

    // make sure we have the ratios with respect to the current and not the initial determinant value
    pcps::xscal(a_ratios.size(), (_init_value/this->_curr_val).get(), &a_ratios[0], 1);
    pcps::xscal(b_ratios.size(), (_init_value/this->_curr_val).get(), &b_ratios[0], 1);
    //pcps::xscal(a_ratios.size(), pcps::unity<S>() / (this->_curr_val).get(), &a_ratios[0], 1);
    //pcps::xscal(b_ratios.size(), pcps::unity<S>() / (this->_curr_val).get(), &b_ratios[0], 1);

  // if we are sampling with the local rdm as the guiding function then we already have the ratios
  } else {
    
    a_ratios = _acti_a_ratios;
    b_ratios = _acti_b_ratios;

  }

}

template void pcps::RestPairingFunc<double>::single_excitation_ratios(const pcps::Configuration<double> & config,
                                                                      std::vector<double> & a_ratios,
                                                                      std::vector<double> & b_ratios);
template void pcps::RestPairingFunc<std::complex<double> >::single_excitation_ratios(
 const pcps::Configuration<std::complex<double> > & config,
 std::vector<std::complex<double> > & a_ratios,
 std::vector<std::complex<double> > & b_ratios);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief Computes the ratios
///        <n| ( p^(+)_a q_a + q^(+)_a p_a ) * ( p^(+)_b q_b + q^(+)_b p_b ) |Psi> / <n|Psi>
///        for the current configuration <n|.
///
/// \param[in]      config      the lattice configuration the object is currently initialized for
/// \param[out]     ratios      on exit, the matrix of pqpq raitios
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S> void pcps::RestPairingFunc<S>::pqpq_ratios(const pcps::Configuration<S> & config, std::vector<S> & ratios) {

  // size the output vector and fill it with zeros
  if (ratios.size() != _nas * _nas)
    ratios.resize(_nas * _nas);
  std::fill(ratios.begin(), ratios.end(), pcps::zero<S>());

  // pre-compute all components of the alpha V (row diff) vectors      slow O(n^2)
  //...

  // pre-compute all components of the beta U (col diff) vectors       slow O(n^2)
  //...

  // pre-compute inverse times all components of the alpha V vectors using one big dgemm call   fast O(n^3)
  //...

  // compute ratios
  for (int i = 0; i < config.a_occ_to_site().size(); i++) {

    const int pa = config.a_occ_to_site().at(i);

    // get the U (delta vec) and UZT vectors related to the alpha p     fast O(n^2)
    //...

    for (int j = 0; j < int(config.a_unocc_to_site().size()); j++) {

      const int qa = config.a_unocc_to_site().at(j);

      // get the beta p and q indices
      int pb, qb;
      if ( config.b_unocc_to_site().at(pa) >= 0 && config.b_occ_to_site().at(qa) >= 0 ) {
        pb = config.b_occ_to_site().at(qa);
        qb = config.b_unocc_to_site().at(pa);
      } else if ( config.b_unocc_to_site().at(qa) >= 0 && config.b_occ_to_site().at(pa) >= 0 ) {
        pb = config.b_occ_to_site().at(pa);
        qb = config.b_unocc_to_site().at(qa);
      } else {
        continue; // skip this alpha pair if there is no matching beta pair
      }

      // build the V (row diff) and VZ vectors related to alpha from their pre-computed components    fast O(n^3)
      //...

      // build the U (col diff) vector related to beta from its pre-computed components    fast O(n^3)
      //...

      // build the V (delta vec) and VZ vectors related to beta    fast O(n^3)
      //...

      // compute the "really easy" inner product deriving from the two delta vecs    fast O(n^2)
      //...

      // set the U and V masks appropriately    fast O(n^2)
      //...

      // compute the determinant ratio    fast O(n^3)
      //...

    }

  }

}

template void pcps::RestPairingFunc<double>::pqpq_ratios(const pcps::Configuration<double> &, std::vector<double> &);
template void pcps::RestPairingFunc<std::complex<double> >::pqpq_ratios(const pcps::Configuration<std::complex<double> > &,
                                                                        std::vector<std::complex<double> > &);

// ///////////////////////////////////////////////////////////////////////////////////////////////////
// /// \brief Exchanges the given rows or columns in the inverse of the occupied pairing matrix.
// ///
// /// \param[in]      b    whether the change is for an alpha or beta update, (0-alpha, 1-beta)
// /// \param[in]      i    the index of the 1st row or column to swap
// /// \param[in]      j    the index of the 2nd row or column to swap
// ///
// ///////////////////////////////////////////////////////////////////////////////////////////////////
// template<class S> void pcps::RestPairingFunc<S>::swap_row_col(const int b, const int i, const int j) {
// 
//   // for alpha, swap columns of the inverse matrix
//   if (b == 0) {
// 
//     for (int p = 0; p < _nap; p++)
//       std::swap(_inv_mat[_nap*p+i], _inv_mat[_nap*p+j]);
// 
//   // for beta, swap rows of the inverse matrix
//   } else {
// 
//     for (int p = 0; p < _nap; p++)
//       std::swap(_inv_mat[_nap*i+p], _inv_mat[_nap*j+p]);
// 
//   }
// 
// }

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief Updates the object so that it is initialized for a new configuration that differs from
///        the currently initialized configuration as specified.
///
/// \param[in]      new_config  the lattice configuration we want to object to be initialized for
/// \param[in]      old_config  the lattice configuration the object is currently initialized for
/// \param[in]      diff        an object specifying how the configurations differ
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S> void pcps::RestPairingFunc<S>::update(const pcps::Configuration<S> & new_config,
                                                        const pcps::Configuration<S> & old_config,
                                                        const pcps::ConfigDiff<S> & diff)
{

  pcps::start_timer("update");

  // check that exactly one particle moved
  if ( diff.newly_empty().size() != 1 || diff.newly_full().size() != 1 )
    throw pcps::Exception("pcps::RestPairingFunc<S>::get_value currently only supports single particle moves");

  const int new_ind = diff.newly_full().at(0);
  const int old_ind = diff.newly_empty().at(0);

  // check that the move did not change the particle's spin
  if ( diff.newly_empty().at(0) % 2 != diff.newly_full().at(0) % 2 )
    throw pcps::Exception("pcps::RestPairingFunc<S>::get_value requires the two sites to be of the same spin");

  // if we are sampling with the wavefunction as the guiding function
  if ( _wfn_sampling ) {

    //this->initialize(new_config);

    this->_curr_val = this->get_value(new_config, old_config, diff);
    _rest_value = this->_curr_val;

    // get the U and V vectors for the single excitation
    {
      _prop_desc = _rest_desc;
      _prop_desc.new_occ_site() = new_ind;
      _prop_desc.old_occ_site() = old_ind;
      this->get_update_vecs(0, _prop_desc);
    }

    if ( new_ind % 2 == 0 ) {
      assert( _Umask.at(0) == 1 );
      assert( _Vmask.at(0) == 0 );
    } else {
      assert( _Umask.at(0) == 0 );
      assert( _Vmask.at(0) == 1 );
    }

    //{
    //  const int rc_ind = _rest_desc.get_row_col_index(old_ind);
    //  _rest_desc.build_occupied_matrix(this->pairing_matrix(), _Zmat);
    //  if ( new_ind % 2 == 0 )
    //    pcps::xaxpy(_nap, pcps::unity<S>(), &_Vmat.at(0), 1, &_Zmat.at(_nap*rc_ind), 1);
    //  else
    //    pcps::xaxpy(_nap, pcps::unity<S>(), &_Umat.at(0), 1, &_Zmat.at(rc_ind), _nap);
    //  S det_val;
    //  pcps::matrix_inverse_lu(_nap, det_val, &_Zmat.at(0), &_ZTmat.at(0), &_iwork.at(0));
    //}

    //{
    //  const int rc_ind = _rest_desc.get_row_col_index(old_ind);
    //  if ( new_ind % 2 == 0 ) {
    //    std::fill(_Umat.begin(), _Umat.end(), pcps::zero<S>());
    //    _Umat.at(rc_ind) = pcps::unity<S>();
    //  } else {
    //    std::fill(_Vmat.begin(), _Vmat.end(), pcps::zero<S>());
    //    _Vmat.at(rc_ind) = pcps::unity<S>();
    //  }
    //}

    // update the inverse matrix using the woodbury matrix identity
    _Imat[0] = pcps::unity<S>();
    _Qmask[0] = 0;
    pcps::woodbury_matrix_identity(_nap, 1, &_Zmat[0], &_Umat[0], &_Vmat[0], &_Imat[0],
                                   &_UZTmat[0], &_VZmat[0], &_Qmat[0],
                                   &_Umask[0], &_Vmask[0], &_Qmask[0], &_work[0], &_iwork[0]);

    // update the descriptors
    _rest_desc.update(new_ind, old_ind);
    _acti_desc = _rest_desc;
    _prop_desc = _rest_desc;

    //_rest_desc.build_occupied_matrix(this->pairing_matrix(), _ZTmat);
    //std::vector<S> eye(_nap*_nap, pcps::zero<S>());
    //pcps::xgemm('N', 'N', _nap, _nap, _nap, pcps::unity<S>(), &_Zmat.at(0), _nap, &_ZTmat.at(0), _nap, pcps::zero<S>(), &eye.at(0), _nap);
    //for (int i = 0; i < _nap; i++)
    //  eye.at(_nap*i+i) -= pcps::unity<S>();
    //const double error = std::sqrt( std::abs( pcps::xdotc(_nap*_nap, &eye.at(0), 1, &eye.at(0), 1) ) );
    //std::cout << boost::format("update error = %20.12f") % error << std::endl;
    //assert( std::abs(error) < 1.0e-6 );

  // if we are sampling with the local rdm as the guiding function
  } else {

//    const int best_p = _best_prop_p;
//    const int best_q = _best_prop_q;
//    pcps::Configuration<S> temp_config(new_config);
//    assert( temp_config.fock_config().at(best_p) == 1 );
//    assert( temp_config.fock_config().at(best_q) == 0 || best_p == best_q );
//    temp_config.fock_config().at(best_p) = 0;
//    temp_config.fock_config().at(best_q) = 1;
//    this->initialize(temp_config);
//    //_rest_desc.new_occ_site() = best_p;
//    //_rest_desc.old_occ_site() = best_q;
//    //_acti_desc = _rest_desc;
//    //_acti_desc.update(_rest_desc.new_occ_site(), _rest_desc.old_occ_site());
//    //_prop_desc = _acti_desc;
//    {
//      pcps::ConfigDiff<S> temp_diff;
//      temp_diff.newly_full().assign(1, best_p);
//      temp_diff.newly_empty().assign(1, best_q);
//      this->get_value(new_config, temp_config, temp_diff);
//      _acti_a_ratios  = _prop_a_ratios;
//      _acti_b_ratios  = _prop_b_ratios;
//      _acti_ratio_sum = _prop_ratio_sum;
//    }
//    _rest_desc.new_occ_site() = best_p;
//    _rest_desc.old_occ_site() = best_q;
//    _acti_desc = _rest_desc;
//    _acti_desc.update(_rest_desc.new_occ_site(), _rest_desc.old_occ_site());
//    _prop_desc = _acti_desc;

    // Here we will use the "best" (largest magnitude of determinant) of the excited
    // configurations surrounding the proposed new active configuration.
    // This new resting configuration should have already been decided upon
    // during the calculation of the local rdm guiding function during the accept/reject step.

    // get the determinant value at the new resting configuration
    _rest_value = _rest_value * _best_prop_ratio;

    // get the ratios and their sum at the new active configuration
    _acti_a_ratios  = _prop_a_ratios;
    _acti_b_ratios  = _prop_b_ratios;
    _acti_ratio_sum = _prop_ratio_sum;

    // prepare the identity matrix for a rank three update
    std::fill(_Imat.begin(), _Imat.end(), pcps::zero<S>());
    _Imat[0*3+0] = pcps::unity<S>();
    _Imat[1*3+1] = pcps::unity<S>();
    _Imat[2*3+2] = pcps::unity<S>();

    // reset the masks
    std::fill(_Umask.begin(), _Umask.end(), 0);
    std::fill(_Vmask.begin(), _Vmask.end(), 0);
    std::fill(_Qmask.begin(), _Qmask.end(), 0);

    // prepare information for the active to proposed move
    assert( _acti_desc.new_occ_site() == new_ind );
    assert( _acti_desc.old_occ_site() == old_ind );

    // get a description of the pairing function at the proposed configuration
    _prop_desc = _acti_desc;
    _prop_desc.update(new_ind, old_ind);
    _prop_desc.new_occ_site() = _best_prop_q;
    _prop_desc.old_occ_site() = _best_prop_p;

    // get the update vectors for the move from resting to active
    this->get_update_vecs(0, _rest_desc);

    // get the update vectors for the move from active to proposed
    this->get_update_vecs(1, _acti_desc);

    // get the update vectors for the move from proposed to the new resting
    this->get_update_vecs(2, _prop_desc);

    // update the inverse matrix for the new resting configuration
    pcps::woodbury_matrix_identity(_nap, 3, &_Zmat[0],
                                   &_Umat[0], &_Vmat[0], &_Imat[0],
                                   &_UZTmat[0], &_VZmat[0], &_Qmat[0],
                                   &_Umask[0], &_Vmask[0], &_Qmask[0],
                                   &_work[0], &_iwork[0]);

    // update the resting configuration descriptor
    _rest_desc.update(_rest_desc.new_occ_site(), _rest_desc.old_occ_site());
    _rest_desc.update(_acti_desc.new_occ_site(), _acti_desc.old_occ_site());
    _rest_desc.update(_prop_desc.new_occ_site(), _prop_desc.old_occ_site());
    _rest_desc.new_occ_site() = _best_prop_p;
    _rest_desc.old_occ_site() = _best_prop_q;

    // update the active configuration descriptor
    _acti_desc.update(_acti_desc.new_occ_site(), _acti_desc.old_occ_site());

    //if ( _rest_desc.new_occ_site() == _rest_desc.old_occ_site() )
    //  this->initialize(new_config);

  }

  // update the transpose of the inverse matrix
  for (int i = 0; i < _nap; i++)
  for (int j = 0; j < _nap; j++)
    _ZTmat[_nap*i+j] = _Zmat[_nap*j+i];

//  // build the little 3x3 matrix for the Woodbury matrix identity
//  this->build_little_mat(diff);
//
//  // update the value of the coefficient
//  this->_curr_val = this->get_single_excite_value(diff.newly_full().at(0), diff.newly_empty().at(0), _little_mat);
//
//  // get whether the move was alpha or beta
//  const int b = diff.newly_empty().at(0) % 2;
//
//  // get the row or column index of the particle that was destroyed
//  const int rc_ind = this->get_row_col_index(diff.newly_empty().at(0));
//
//  // get the row or column difference
//  this->get_row_col_diff(diff, &_work_mat.at(0));
//
//  // construct the big ( inverse times U ) and ( V times inverse ) matrices
//  std::fill(_bigUI.begin(), _bigUI.end(), pcps::zero<S>());
//  std::fill(_bigVI.begin(), _bigVI.end(), pcps::zero<S>());
//  if ( b == 0 ) {
//
////    S x = pcps::unity<S>();
////    for (int i = 0; i < _nap; i++)
////      x += _work_mat.at(i) * _inv_mat.at(_nap*i+rc_ind);
////    for (int i = 0; i < _nap; i++) {
////      _more_work.at(i) = pcps::zero<S>();
////      for (int j = 0; j < _nap; j++)
////        _more_work.at(i) += _work_mat.at(j) * _inv_mat.at(_nap*j+i);
////    }
////    for (int i = 0; i < _nap; i++)
////    for (int j = 0; j < _nap; j++)
////      _work_mat.at(_nap*i+j) = _inv_mat.at(_nap*i+rc_ind) * _more_work.at(j);
////    for (int i = 0; i < _nap; i++)
////    for (int j = 0; j < _nap; j++)
////      _inv_mat.at(_nap*i+j) -= _work_mat.at(_nap*i+j) / x;
//
//    assert( diff.noise_change_codes().at(0) != 'u');
//    assert( diff.noise_change_codes().at(1) != 'u');
//
//    // UI
//    pcps::xcopy(_nap, &_inv_mat.at(rc_ind), _nap, &_bigUI.at(_nap*0), 1);
//    if ( diff.noise_change_codes().at(0) == 'v' )
//      pcps::xcopy(_nap, &_Iu1.at(0), 1, &_bigUI.at(_nap*1), 1);
//    if ( diff.noise_change_codes().at(1) == 'v' )
//      pcps::xcopy(_nap, &_Iu2.at(0), 1, &_bigUI.at(_nap*2), 1);
//
//    // VI
//    pcps::xgemm('N', 'T', 1, _nap, _nap, pcps::unity<S>(), &_work_mat.at(0), 1,
//                &_inv_mat.at(0), _nap, pcps::zero<S>(), &_bigVI.at(_nap*0), 1);
//    if ( diff.noise_change_codes().at(0) == 'v' )
//      pcps::xgemm('N', 'T', 1, _nap, _nap, pcps::unity<S>(), &diff.noise_changes().at(0).at(0), 1,
//                  &_inv_mat.at(0), _nap, pcps::zero<S>(), &_bigVI.at(_nap*1), 1);
//    if ( diff.noise_change_codes().at(1) == 'v' )
//      pcps::xgemm('N', 'T', 1, _nap, _nap, pcps::unity<S>(), &diff.noise_changes().at(1).at(0), 1,
//                  &_inv_mat.at(0), _nap, pcps::zero<S>(), &_bigVI.at(_nap*2), 1);
//    
//  } else {
//
////    S x = pcps::unity<S>();
////    for (int i = 0; i < _nap; i++)
////      x += _work_mat.at(i) * _inv_mat.at(_nap*rc_ind+i);
////    for (int i = 0; i < _nap; i++) {
////      _more_work.at(i) = pcps::zero<S>();
////      for (int j = 0; j < _nap; j++)
////        _more_work.at(i) += _work_mat.at(j) * _inv_mat.at(_nap*i+j);
////    }
////    for (int i = 0; i < _nap; i++)
////    for (int j = 0; j < _nap; j++)
////      _work_mat.at(_nap*i+j) = _inv_mat.at(_nap*rc_ind+j) * _more_work.at(i);
////    for (int i = 0; i < _nap; i++)
////    for (int j = 0; j < _nap; j++)
////      _inv_mat.at(_nap*i+j) -= _work_mat.at(_nap*i+j) / x;
//
//    assert( diff.noise_change_codes().at(0) != 'v');
//    assert( diff.noise_change_codes().at(1) != 'v');
//
//    // UI
//    pcps::xgemm('T', 'N', _nap, 1, _nap, pcps::unity<S>(), &_inv_mat.at(0), _nap,
//                 &_work_mat.at(0), _nap, pcps::zero<S>(), &_bigUI.at(_nap*0), _nap);
//    if ( diff.noise_change_codes().at(0) == 'u' )
//      pcps::xgemm('T', 'N', _nap, 1, _nap, pcps::unity<S>(), &_inv_mat.at(0), _nap,
//                  &diff.noise_changes().at(0).at(0), _nap, pcps::zero<S>(), &_bigUI.at(_nap*1), _nap);
//    if ( diff.noise_change_codes().at(1) == 'u' )
//      pcps::xgemm('T', 'N', _nap, 1, _nap, pcps::unity<S>(), &_inv_mat.at(0), _nap,
//                  &diff.noise_changes().at(1).at(0), _nap, pcps::zero<S>(), &_bigUI.at(_nap*2), _nap);
//
//    // VI
//    pcps::xcopy(_nap, &_inv_mat.at(_nap*rc_ind), 1, &_bigVI.at(_nap*0), 1);
//    if ( diff.noise_change_codes().at(0) == 'u' )
//      pcps::xcopy(_nap, &_v1I.at(0), 1, &_bigVI.at(_nap*1), 1);
//    if ( diff.noise_change_codes().at(1) == 'u' )
//      pcps::xcopy(_nap, &_v2I.at(0), 1, &_bigVI.at(_nap*2), 1);
//
//  }
//
//  // invert the little matrix
//  {
//    S det_val;
//    pcps::matrix_inverse_lu(3, det_val, &_little_mat.at(0), &_more_work.at(0), &_iwork.at(0));
//  }
//
//  // contract the little matrix's inverse with the big VI matrix
//  pcps::xgemm('N', 'N', _nap, 3, 3, pcps::unity<S>(), &_bigVI.at(0), _nap,
//              &_little_mat.at(0), 3, pcps::zero<S>(), &_more_work.at(0), _nap);
//
//  // contract the intermediate with the big UI matrix to update the inverse matrix
//  pcps::xgemm('N', 'T', _nap, _nap, 3, -pcps::unity<S>(), &_more_work.at(0), _nap,
//              &_bigUI.at(0), _nap, pcps::unity<S>(), &_inv_mat.at(0), _nap);
//
//  //// reorder rows or columns so that the occupied site indices are in ascending order again
//  //std::vector<int> & occ_sites = ( b == 0 ? _a_sites : _b_sites );
//  //occ_sites.at(rc_ind) = diff.newly_full().at(0) / 2;
//  //for (int i = rc_ind + 1; i < _nap && occ_sites[i] < occ_sites[i-1]; i++) {
//  //  std::swap(occ_sites[i], occ_sites[i-1]);
//  //  this->swap_row_col(b, i-1, i);
//  //}
//  //for (int i = rc_ind - 1; i >= 0 && occ_sites[i] > occ_sites[i+1]; i--) {
//  //  std::swap(occ_sites[i], occ_sites[i+1]);
//  //  this->swap_row_col(b, i+1, i);
//  //}
//
//  // update the occupancy list
//  std::vector<int> & occ_sites = ( b == 0 ? _a_sites : _b_sites );
//  assert( occ_sites.at(rc_ind) == diff.newly_empty().at(0) / 2 );
//  occ_sites.at(rc_ind) = diff.newly_full().at(0) / 2;
//
//  // update the permutation that keeps track of the ordering of the rows' or columns' sites
//  if ( b == 0 )
//    this->initialize_permutation(_a_sites, _a_ordered_perm, _a_perm_pos);
//  else
//    this->initialize_permutation(_b_sites, _b_ordered_perm, _b_perm_pos);
//
//  // compute the products of the inverse times the noise vectors
//  this->compute_noise_products(new_config);
//
//  //this->print_occupancy_info();
//
//  if (false) {
//    std::cout << "checking noise removal after update:" << std::endl << std::endl;
//    this->test_noise_removal(new_config);
//  }

  pcps::stop_timer("update");

}

template void pcps::RestPairingFunc<double>::update(const pcps::Configuration<double> &,
                                                    const pcps::Configuration<double> &,
                                                    const pcps::ConfigDiff<double> &);
template void pcps::RestPairingFunc<std::complex<double> >::update(const pcps::Configuration<std::complex<double> > &,
                                                                   const pcps::Configuration<std::complex<double> > &,
                                                                   const pcps::ConfigDiff<std::complex<double> > &);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  retrieves the correct guiding function value for the n2 sto3g testing
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S> S pcps::RestPairingFunc<S>::correct_value(const pcps::Configuration<S> & config, const char val_type) const {

  std::stringstream key_stream;
  for (int i = 0; i < config.fock_config().size(); i+=2)
    key_stream << config.fock_config().at(i);
  for (int i = 1; i < config.fock_config().size(); i+=2)
    key_stream << config.fock_config().at(i);
  const std::string key_str = key_stream.str();

  std::map<std::string, std::pair<double, double> >::const_iterator it = _cmap.find(key_str);

  if ( it == _cmap.end() )
    throw pcps::Exception("could not find correct value for configuration %s") % key_str;

  //if ( !_wfn_sampling )
  //  throw pcps::Exception("GOTCHA");

  assert( val_type == 'v' || val_type == 'e' );

  return pcps::unity<S>() * ( val_type == 'v' ? it->second.first : it->second.second );

}

template double pcps::RestPairingFunc<double>::correct_value(const pcps::Configuration<double> &, const char) const;
template std::complex<double>  pcps::RestPairingFunc<std::complex<double> >::correct_value(
 const pcps::Configuration<std::complex<double> > &, const char) const;

